Understanding how spatial variation is linked to diversity maintenance in natural communities is a pillar of plant community ecology. Theoretically, a variable landscape can maintain diversity via niche partitioning: different species can trade off in performing better or worse depending on the conditions of the patch they are growing in, and as a result, more species can sustainably coexist in a community than if it were spatially heterogeneous. In the hyperdiverse system of native annual plants in Western Australia, fallen logs may be one of the greatest contributors to generating spatial variation that could help maintain species diversity. Considerable anecdotal evidence suggests that fallen logs generate spatial variation, or patchiness, in the environment (Figure 1), and that species or assemblages of plants may respond differently depending on if they are near logs or not. Despite such anecdotal evidence, it is yet unknown if and how fallen logs contribute to maintaining species diversity in the native annual plant communities of the Western Australian wheat belt.

Figure 1: image of annual plant halos around logs
Figure 1: image of annual plant halos around logs

The project will address the following questions:

Q1) Are/how are plant communities in fallen log patches different from patches that are in the open?

Q2) Why are plant communities in fallen log patches different from patches in the open?

Q3) Are/how are plant species performances affected by proximity to fallen logs?

Hypotheses

The null hypothesis, H0, is that annual plants in fallen log patches are not different in diversity, abundance, or composition from open patches.

In addition to the null hypothesis, the following constitute four, non-mutually exclusive hypotheses concerning how fallen logs may introduce spatial variation in the environment. I include corresponding predictions for how plant communities may differ between fallen log patches as compared to open patches.

H1: Log decomposition creates islands of fertility directly around the fallen log.
Prediction 1: Nutrient composition around logs will be higher than in open plots

Prediction 2: Variations in nutrient composition in log vs open environments will correspond to variations in species composition, abundance, and/or richness in these environments.

Prediction 3: All sown plants will perform best in environments where organic logs have been left ‘insitu’. In locations where logs have been removed or replaced with pvc, the legacy of the nutrient island effect will yield higher sown plant performance than when compared to locations where logs have never been. The effect of the nutrient island in locations where logs have been added to open environments should yeild higher plant performance over time. note: performance is measured in terms of germination rate, survival to fruiting, fecundity, and/or biomass.

H2: Fallen logs alter the microclimate directly around them by providing shade.
Prediction 1: Shade and temperature around logs vs in open plots will be different

Prediction 2: Variation in shade and temperature in log vs open environments will correspond to variation in species composition, abundance, and/or richness in these environments.

Prediction 3: All sown plants will perform best in environments where there are organic or pvc logs, no matter if they have been recently moved or not.

H3: Fallen logs trap dispersing seeds as they are blown along the ground.

Prediction 1: Dispersing seeds accumulate around logs, leading to a denser stand of plants in fallen log patches. Plant abundance in fallen log patches will be higher as compared to open patches. Rare plants will be more common in fallen log patches as compared to open patches

Prediction 2: All sown plants will perform the same in all experimental environments

H4: At least some species perform differently according to variation in log vs. open environments and have short dispersal kernels, causing fitness-density covariance

Figure 2: Photo before germination, after a rain. Notice the seeming wet halo under and around the branch
Figure 2: Photo before germination, after a rain. Notice the seeming wet halo under and around the branch


Experimental Design

In this experiment, 224 plots are arranged in 7 blocks of 32 plots each within the Caron Dam nature reserve. A map can be found here.
note: the location info for 3.02 is probably incorrect as of May 2022, and location info is currently unavailable for plots 6.25 and 7.19

Each block is approximately 30m X 30m in area. Plots are 1m long and linear, and have a pin tag on either end (see Figure 3). The pin tags have the identity of the plot written on them in the form of “blocknumber.plotnumber”. Plots are 1m or more away from each other.

In each block, plot environments can be one of six types:
- A 1m log that is out in the open (open_with_log, 4 plots)
- A 1m log that is a part of a tree (insitu_log, 4 plots)
- A 1m pvc pipe that is out in the open (open_with_pvc, 4 plots)
- A 1m pvc pipe that is a part of a tree (insitu_pvc, 4 plots)
- A plot that is out in the open (open, 8 plots)
- A gap in a log where a log used to be (gap, 8 plots)

In half of the plots (not including open plots), the addition, exchange, or removal of logs or pvc to the environment was implemented in October 2020, before seed dispersal. In the other half of these plots, these manipulations were implemented after seed dispersal, in March 2021.

Within each 1m long plot, there is a ~20cm long microtransect. The ends of the microtransects are marked by a nail and a washer sunken into the ground. Each microtransect is approximately 21 cm in internal length from inner washer edge to inner washer edge. Microtransects are not sided.

In half of all plots, seeds were sown in March 2021 and February 2022. In these plots, 15 seeds each of Trachymene ornata (TROR), Goodenia rosea (GORO), and Trachymene cyanopetala (TRCY) are sown outside of the microtransects as in the diagram. These plants were selected because they represent plants common to communities next to logs (TROR), out in the open (GORO), or both (TRCY). The plots where seeds were sown are called ‘lambda’ plots as noted in Figure 3. In the dataset, the rows with a “1” in the ‘seeding_trt’ column are the plots that had seeds sown into them.

Figure 3: plot schematic
Figure 3: plot schematic

Datasets

The sets of data that we have collected for this arm of the project are the following.

  1. Community data, before and after the experiment was implemented.
  1. Soil nutrient analysis in the open and insitu log plots.
  1. Performance data of TROR, TRCY, and TROR.

Methods and analysis

Q1 Are/how are plant communities in fallen log patches different from patches that are in the open?

Overview of results

Quick recap:

  1. Plant abundance does not differ between log and open plots across all years.

  2. Plant diversity is higher in log patches in 2022 but no significant difference was detected in 2020 and 2021.

  3. Plot type (log vs open) explains a near-zero portion of the variance in plant species composition across plots across years.

To compare the characteristic of plant communities in fallen log patches and open patches, we assessed differences in plant abundance, diversity using the Shannon diversity index, and plant composition summarized at the transect level. The analysis is organized by year - we analysed the plot-level community data for each year separately: 2020, 2021, and 2022. In other words, each row in the data set represents an individual plot from a specific year. In these analyses, we evaluated three response variables: abundance, diversity, and plant composition.

2020

There are three response variables: count (Poisson), Shannon diversity index (zero-inflated glmer with Gaussian), and composition matrix.

  1. Abundance – No difference in total plant count between log patches and open patches (p = 0.355). Model: total count ~ treatment(log/open) + (1|block).

  2. Diversity – Plant diversity (SDI) is marginally lower in open patches than log patches (p = 0.063). Model: diversity ~ treatment(log/open) + (1|block) with zero-inflation terms=~1

  3. Composition – Partial rda partitions variance in plant composition by R-squared. As expressed by the constrained term, ‘treatment’ (p = .001) explains only 1.36% of variations in plant composition. The conditional term ‘block’ explains 9.17% of variance in plant composition. 89.5% of the variance in plant composition remains unexplained in unconstrained terms. The analysis of variation partitioning uses an adjusted R-squared method to express the variance explained by the explanatory matrix. The variable ‘treatment’ alone explains about 1% of the total variance, while the ‘block’ explains roughly 7% of the variance. Note that the shared fraction in the Venn diagram is not shown as it is negative, likely due to ‘treatment’ being a suppressor variable with a close-to-zero relationship to plant composition. (Read: https://www.jstor.org/stable/20069271). In a parallel set of analyses, the 4-D NMDS ordination has stress values of 4.8%, suggesting a good fit of ordination. ANOSIM indicates no difference in plant composition between open and log plots (P = 0.5985).

2021

There are three response variables: count (Poisson), Shannon diversity index (zero-inflated glmer with Gaussian), and composition matrix.

  1. Abundance – No difference between total plant count between log patches and open patches (p = 0.322). Model: total count ~ treatment(log/open) + (1|block).

  2. Diversity – No difference in plant diversity (SDI) between open patches and log patches (p = 0.822). Note that the zero-inflated model has a significant constant term. Model: diversity ~ treatment(log/open) + (1|block) with zero-inflation terms=~1

  3. Composition – As expressed by the constrained term, ‘treatment’ (p = .01) explains only 2.15% of variations in plant composition. The conditional term ‘block’ explains 13.5% of variance in plant composition. 84.4% of the variance in plant composition remains unexplained in unconstrained terms. In terms of adjusted R-squared, the Venn diagram shows that ‘treatment’ alone explains about 1% and ‘block’ explains roughly 7% of the total variance in plant composition. Similarly, ‘treatment’ is likely to be a suppressor. In a parallel set of analyses, the 5-D NMDS ordination has stress values of 4.1%, suggesting a very good fit of ordination. ANOSIM indicates a significant difference in plant composition between open and log plots (P = 0.0088) but the magnitude of difference is rather small (R = 0.28 - ratio of dissimilarity in plant composition between log and open plot types).

2022

There are three response variables: count (Poisson), Shannon diversity index (zero-inflated glmer with Gaussian), and composition matrix.

  1. Abundance – No difference between total plant count between log patches and open patches (p = 0.322). Model: total count ~ treatment(log/open) + (1|block).

  2. Diversity – Plant diversity (SDI) is significantly lower in open patches than log patches (p < 0.001). Model: diversity ~ treatment(log/open) + (1|block) with zero-inflation terms=~1

  3. Composition – As expressed by the constrained term, ‘treatment’ (p = 0.002) explains only 3.08% of variations in plant composition. The conditional term ‘block’ explains 12.7% of variance in plant composition. 84.2% of the variance in plant composition remains unexplained in unconstrained terms. In terms of adjusted R-squared, the Venn diagram shows that ‘treatment’ alone explains about 2% and ‘block’ explains roughly 6% of the total variance in plant composition. Similarly, ‘treatment’ is likely to be a suppressor. In a parallel set of analyses, the 4-D NMDS ordination has stress values of 4%, suggesting a very good fit of ordination. ANOSIM indicates a significant difference in plant composition between open and log plots (P = 0.0017) but the magnitude of difference is rather small (R = 0.31 - ratio of dissimilarity in plant composition between log and open plot types).

Statistical Methods

  1. Abundance comparison

Comparison of plant abundance between log and open plots involves two types of data: (1) the total plant counts of individual plots of all plots in 2020 before the experiment set up, which includes two initial treatments - “log” and “open”; and (2) the total plant counts from insitu_log and insitu_open plots in 2021 and 2022 respectively. All count data used for abundance comparison includes species with unknown identity. We used a generalized linear mixed-effects model (GLMER) to model the abundance of plant individuals (2020, 2021 and 2022 data) in unaltered fallen log patches and open patches using a Poisson distribution structure with ‘block’ as a random term.

  1. Diversity comparison

A 20-centimeter linear transect was set up for each plot in 2020. We recorded species identity and individual counts for all species occurring along the transects in 2020, 2021 and 2022 respectively. Comparison of plant diversity between log and open patches involves two types of data: 1. The count of plants for each species in individual plots of all plots in the year 2020, before the experiment setup. The data is differentiated based on the initial treatment of “log” and “open”. 2. The count of plants for each species in individual plots from insitu_log and insitu_open plots in the years 2021 and 2022, respectively. All species with unknown identity have been excluded from the data used for diversity comparison.

The diversity of plants in each community (i.e. plot) was measured by Shannon diversity index (SDI). The distribution of SDI values is zero-inflated as testing with DHARMa. Note that a zero SDI could mean 0 plant in a community or one species without any neighbours. To account for the zero inflation, we used glmmTMB to model the diversity of plant community. We assumed that zeros are generated by different mechanisms in our model and the probability of having an extra zero is uniform for all observations by setting zi=~1.

  1. Composition

Comparison of plant composition between log patches and open patches involves two types of data: 1. The count of plants for each species in individual plots of all plots in the year 2020, before the experiment setup. The data is differentiated based on the initial treatment of “log” and “open”. 2. The count of plants for each species in individual plots from insitu_log and insitu_open plots in the years 2021 and 2022, respectively. All data used for diversity comparison excludes species with unknown identity. For transects where no plants were found, an artificial species, “x,” was added to the species composition matrix to represent zero plants. Data is structured at the plot level with each row representing the composition of plant species in an individual plot. There are 223 plots in 2020, 86 plots in 2021, and 84 plots in 2022, evenly nested in 7 blocks.

We compared plant species composition between plot treatments (i.e. open vs log) with two separate streams of analyses (1) partial Redundancy Analysis (RDA) + permutation test (ANOVA.acc)+ variation partioning analysis (Venn diagram) and (2) NMDS + ANOSIM (Analysis of Similarities).

(3.1) We ran a partial Redundancy Analysis (RDA) to test the correlation between species composition (squared-root transformed) and plot treatment (log vs open patch) with ‘block’ as a conditional term. The effect of plot treatment on plant species composition was adjusted for the covariate ‘block’. Significance of model and ‘treatment’ term were tested with a permutation test (anova.cca function from vegan package). Variation in plant composition explained by plot treatment and block were partitioned and reported in an adjusted R-squared value in a Venn diagram.

(3.2) Alternatively, we used Non-Metric Dimensional Scaling (NMDS) to analyse the differences in plant communities with Bray-Curtis dissimilarity metrics. For each combination of NMDS axes generated, we plotted a two-dimensional NMDS ordination plot to visualize and interpret the relative position of plant community between plot types and blocks. To better understand the contribution of each species to the compositional changes, we extracted and plotted the ordination scores for each species along each respective NMDS axis. These scores represent the weighted average of a species’ abundance score in a sample community along each selected NMDS axes. To further test the degree of similarity in plant species composition between plot treatment types, we used the ANOSIM function, and the resulting R-statistic value depicts the proportion of dissimilarity.

Analysis

Abundance analysis
** 2020 - 2022 **
  Abundance:2020-2022 Abundance:2020 Abundance:2021 Abundance:2022
Predictors Incidence Rate Ratios CI p Incidence Rate Ratios CI p Incidence Rate Ratios CI p Incidence Rate Ratios CI p
(Intercept) 6.79 4.88 – 9.46 <0.001 4.74 4.34 – 5.18 <0.001 7.14 6.20 – 8.22 <0.001 8.82 7.56 – 10.30 <0.001
init [open] 1.00 0.93 – 1.09 0.913 0.94 0.83 – 1.07 0.355 1.09 0.92 – 1.28 0.322 1.04 0.89 – 1.21 0.642
Random Effects
σ2 0.14 0.20 0.12 0.10
τ00 0.00 block 0.00 block 0.00 block 0.01 block
0.08 time      
ICC 0.38 0.00 0.03 0.12
N 7 block 7 block 7 block 7 block
3 time      
Observations 391 221 86 84
Marginal R2 / Conditional R2 0.000 / 0.378 0.004 / 0.009 0.012 / 0.037 0.002 / 0.127

Diversity analysis
Shannon diversity index (SDI): 2020
  SDI [2020]
Predictors Estimates CI p
Count Model
(Intercept) 0.78 0.68 – 0.88 <0.001
init [open] -0.13 -0.26 – 0.01 0.063
(Intercept) 0.51 0.46 – 0.56
Zero-Inflated Model
(Intercept) -20.98 -13594.28 – 13552.31 0.998
Random Effects
σ2 0.51
τ00 block 0.00
ICC 0.01
N block 7
Observations 223
Marginal R2 / Conditional R2 0.008 / 0.013
Shannon diversity index (SDI): 2021
  SDI [2021]
Predictors Estimates CI p
Count Model
(Intercept) 1.06 0.92 – 1.19 <0.001
init [open] -0.02 -0.20 – 0.16 0.822
(Intercept) 0.36 0.28 – 0.45
Zero-Inflated Model
(Intercept) -2.47 -3.46 – -1.48 <0.001
Random Effects
σ2 0.36
τ00 block 0.00
N block 7
Observations 86
Marginal R2 / Conditional R2 0.000 / NA
Shannon diversity index (SDI): 2022
  SDI [2022]
Predictors Estimates CI p
Count Model
(Intercept) 1.21 1.03 – 1.40 <0.001
init [open] -0.38 -0.57 – -0.18 <0.001
(Intercept) 0.42 0.36 – 0.50
Zero-Inflated Model
(Intercept) -22.39 -24824.98 – 24780.19 0.999
Random Effects
σ2 0.42
τ00 block 0.02
ICC 0.04
N block 7
Observations 84
Marginal R2 / Conditional R2 0.068 / 0.105

Composition analysis

Composition dissimilarity 2020

Partitioning_of_variance Inertia Proportion Note Permutation_test Df Variance F P
rda(formula = ass.rel.t0 ~ init + Condition(block)) Model signif. = ***
Total 0.86477 1 adj.r.squared = 0.009696575 treatment 1 0.01176 3.2671 0.001 ***
Conditioned 0.07929 0.09169 Variance explained by ‘block’ Residual 215 0.77372
Constrained 0.01176 0.0136 Variance explained by ‘treatment’
Unconstrained 0.77372 0.89471 Unexplained variance 

## 
## Call:
## anosim(x = assemblies_t0, grouping = blocksum$init, permutations = 9999,      distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: -0.0311 
##       Significance: 0.5919 
## 
## Permutation: free
## Number of permutations: 9999

Composition dissimilarity 2021

Partitioning_of_variance Inertia Proportion Note Permutation_test Df Variance F P
rda(formula = ass.rel.t1 ~ init + Condition(block)) Model signif. = **
Total 0.82992 1 adj.r.squared = 0.01150461 treatment 1 0.01785 1.9885 0.01 **
Conditioned 0.11182 0.13473 Variance explained by ‘block’ Residual 78 0.70025
Constrained 0.01785 0.02151 Variance explained by ‘treatment’
Unconstrained 0.70025 0.84376 Unexplained variance 

## 
## Call:
## anosim(x = assemblies_t1, grouping = blocksum$init, permutations = 9999,      distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.2809 
##       Significance: 0.008 
## 
## Permutation: free
## Number of permutations: 9999

Composition dissimilarity 2022

Partitioning_of_variance Inertia Proportion Note Permutation_test Df Variance F P
rda(formula = ass.rel.t2 ~ init + Condition(block)) Model signif. = **
Total 0.8167 1 adj.r.squared = 0.02119015 treatment 1 0.02511 2.7732 0.002 **
Conditioned 0.10345 0.12667 Variance explained by ‘block’ Residual 76 0.68814
Constrained 0.02511 0.03075 Variance explained by ‘treatment’
Unconstrained 0.68814 0.84259 Unexplained variance 

## 
## Call:
## anosim(x = assemblies_t2, grouping = blocksum$init, permutations = 9999,      distance = "bray") 
## Dissimilarity: bray 
## 
## ANOSIM statistic R: 0.3105 
##       Significance: 0.0012 
## 
## Permutation: free
## Number of permutations: 9999

Q2 Why are plant communities in fallen log patches different from patches in the open?

Overview of results

  1. Soil K and P have similar PCA loadings. In subsequent analysis, we may simplify the model structure by keeping only either K or P and still be able to infer the effects of the other.

  2. Total soil C (Wilcoxon rank sum test, p = 0.025) and plant available N (Wilcoxon rank sum test, p = 0.031) are significantly higher in log patches than open patches

  3. Plant abundance and diversity have a variable relationship with different nutrient elements in different years.

  4. No significant relationships have been found between the composition of plants and soil nutrient tested by post-hoc permutation analysis on RDA. However, it is worth noting that the data may not represent the true effect of soil nutrients on plant composition, as the degree of freedom for each soil nutrient observation is quite small (df=1). Therefore, it is difficult to disentangle the effects of each soil nutrient element on the overall plant species composition. On the other hand, soil nutrient variables are plotted on the NMDS plots by regressing each variable to the unconstrained axes. It means that we cannot use these projections on the NMDS plots to infer the explanatory effects of nutrient elements on plant species composition.

2020

There are three response variables: count (Poisson), Shannon diversity index (zero-inflated glmer with Gaussian), and composition matrix.

  1. Abundance – The best model suggested by model dredging includes all nutrient elements except soil K as additive fixed factors. According to this inclusive model, the abundance of plants is significantly affected by soil N (p=0.002), P (p=0.037), C (p=0.036), and CEC (p=0.020). However, the direction of slopes of each element varies.

  2. Diversity – The best model suggested by model dredging includes base cation (CEC), K and N as additive fixed factors. The Shannon diversity index in the plant community tends to be significantly higher with higher levels of plant-available nitrogen (p=0.001) and potassium (p=0.003).

  3. Composition – Permutation test performed on Redundancy analysis (RDA) shows that there is no significant correlation between any of the nutrient elements and the composition of plant species.

2021

There are three response variables: count (Poisson), Shannon diversity index (zero-inflated glmer with Gaussian), and composition matrix.

  1. Abundance – The best model suggested by model dredging includes only K, N and P as additive fixed factors. According to this inclusive model, the abundance of plants is only significantly affected by soil P (p=0.009) suggesting that plant abundance tends to increase in phosphorus rich soils.

  2. Diversity – The best model suggested by model dredging includes base cation (CEC) and P as additive fixed factors. The Shannon diversity index in the plant community tends to be significantly higher with higher levels of potassium (p=0.006).

  3. Composition – Permutation test performed on Redundancy analysis (RDA) shows that there is no significant correlation between any of the nutrient elements and the composition of plant species. There is only a very weak explanatory power observed in the base cation (CEC: p=0.099)

2022

There are three response variables: count (Poisson), Shannon diversity index (zero-inflated glmer with Gaussian), and composition matrix.

  1. Abundance – The best model suggested by model dredging includes only P and C as additive fixed factors. According to this inclusive model, the abundance of plants is only significantly affected by soil C (p=0.034) suggesting that plant abundance tends to increase in carbon rich soils.

  2. Diversity – The best model suggested by model dredging includes N, K and P as additive fixed factors. The Shannon diversity index in the plant community tends to be significantly higher with higher levels of plant-available nitrogen (p=0.026) and potassium (p=0.004) but lower levels of soil P (p=0.009).

  3. Composition – Permutation test performed on Redundancy analysis (RDA) shows that there is no significant correlation between any of the nutrient elements and the composition of plant species.

Statistical Methods

  1. We performed a principal component analysis using a correlation matrix to examine the multivariate relationships among six different soil properties. These soil properties are soil total carbon (C), potassium (K), phosphorus (P), plant available nitrogen (N), base cation (CEC), and pH.

  2. We used the Wilcoxon rank sum test to compare soil nutrient levels between log and open plots.

  3. Effect of soil nutrient on plant abundance.

This analysis involves three types of data: (1) plant counts of all plots in 2020 (N = 221) before the experiment setup differentiated by the initial treatment “log” and “open”; (2) plant counts from insitu_log and insitu_open plots in 2021 (N = 86) and 2022 (N = 84) respectively, including species with unknown identity; and (3) soil nutrient data collected from one random log plot and one random open plot in each block (N = 14).

Since independent nutrient variables were sampled at block level, this part of the analysis does not include block as the random term to avoid singularity issues. Nutrient data is paired to abundance data summarize at plot level in this analysis. Although the estimate of the slope may remain unaffected, the strength of correlation (i.e. p value) may be magnified with an increase in N, which is a calculation artefact. One of the proposed solutions is to recalculate plant abundance at block level to improve the presentation of data. However, we should pay attention to the df (i.e. df =1) of each nutrient element and decide whether the result is representative of the effect (We should have a discussion on this).

We used a generalized linear model (GLM) with a Poisson distribution structure to model the effect of all soil components on the abundance of plant individuals (data from 2020, 2021, and 2022 respectively). Additionally, we performed a model dredging (using the MuMln package) to select the best model (delta =0) with a subset of soil variables that best describes the distribution of plant abundance. We approach the problem without bias by not holding prior assumptions about the important soil nutrient of interest.

  1. Effect of soil nutrient on plant diversity (Shannon diversity index = SDI)

This analysis involves three types of data: (1) plot-level plant count for each species of all plots in 2020 before the experiment setup differentiated by the initial treatment “log” and “open”; and (2) plot-level plant count for each species from insitu_log and insitu_open plots in 2021 and 2022 respectively. All data used for diversity comparison excludes species with unknown identity. (3) soil nutrient data collected from one random log plot and one random open plot in each block (N = 14). (*Similar issue as described in 3.)

We used the Shannon diversity index to describe the plant species diversity for each transect at the plot level. Similarly to the approach used in Q1, we used a zero-inflated glmer with a Gaussian distribution function to model the effect of all soil properties as additive terms on plant diversity. We did not include ‘block’ as a random term in the model, but the function glmmTMB is able to handle models without a random term. Similarly, we performed a model dredging (using the MuMln package) to select the best model (delta =0) with a subset of soil variables that best describes the distribution of plant diversity.

  1. Effect of soil nutrient on plant composition

We performed NMDS analyses similar to those described in the previous section. However, to enhance the visualization of data, we decided to present the plant species composition data summarized at the block level (N=14). This means that we will have one value for log and one value for open within each block. Additionally, this approach will offer an alternative perspective on the issue highlighted in points 3 and 4. We utilized envfit analyses from the Vegan package to fit soil components onto the nondimensional space. However, it is important to note that soil variables are regressed onto the unconstrained axes in the nondimensional space, and this does not necessarily represent their explanatory power.

Similarly, as an alternative approach, we ran a Redundancy Analysis (RDA) to test the correlation between species composition (squared-root transformed) and all soil variables as additive terms. To avoid pseudo-replication, we excluded ‘block’ as either fixed effect or random effect. We tested the significance of the model and each soil variable term by using a permutation test (anova.cca function from vegan package). We plotted the first two RDA axes to visualize the impact of soil variables on plant species composition.

Analysis

PCA

Nutrient composition comparison between log and open

##   element    w    p_value
## 1       N 41.5 0.03126054
## 2       P 26.5 0.82963146
## 3       K 26.0 0.90151515
## 4       C 42.5 0.02518656
## 5      pH 29.5 0.56096176
## 6     CEC 30.0 0.53496503
Abundance ~ nutrient composition
Abundance ~ soil nutrient [2022]
  Best model Selected inclusive model (2022)
Predictors Log-Mean std. Error CI p Log-Mean std. Error CI p
(Intercept) 2.18 0.04 2.11 – 2.26 <0.001 2.18 0.04 2.10 – 2.26 <0.001
N -0.15 0.05 -0.24 – -0.06 0.002 -0.18 0.06 -0.30 – -0.07 0.002
P -0.09 0.04 -0.17 – -0.01 0.037 -0.07 0.05 -0.16 – 0.02 0.149
C 0.14 0.07 0.01 – 0.26 0.036 0.16 0.07 0.02 – 0.29 0.021
CEC -0.14 0.06 -0.25 – -0.02 0.020 -0.13 0.06 -0.25 – -0.01 0.032
K -0.05 0.06 -0.17 – 0.06 0.336
Observations 84 84
R2 Nagelkerke 0.236 0.246
Abundance ~ soil nutrient [2021]
  Best model Inclusive model
Predictors Log-Mean std. Error CI p Log-Mean std. Error CI p
(Intercept) 2.00 0.04 1.91 – 2.08 <0.001 2.00 0.04 1.92 – 2.08 <0.001
K -0.08 0.05 -0.18 – 0.02 0.100 -0.09 0.06 -0.22 – 0.03 0.154
N -0.09 0.05 -0.19 – 0.01 0.073 -0.10 0.06 -0.23 – 0.02 0.098
P 0.11 0.04 0.03 – 0.19 0.009 0.12 0.04 0.03 – 0.20 0.009
C 0.04 0.08 -0.12 – 0.19 0.598
CEC -0.03 0.07 -0.16 – 0.10 0.664
Observations 86 86
R2 Nagelkerke 0.128 0.132
Abundance ~ soil nutrient [2020]
  Best model Inclusive model
Predictors Log-Mean std. Error CI p Log-Mean std. Error CI p
(Intercept) 1.52 0.03 1.46 – 1.58 <0.001 1.52 0.03 1.46 – 1.58 <0.001
P 0.06 0.03 -0.00 – 0.12 0.054 0.06 0.04 -0.02 – 0.13 0.153
C 0.07 0.03 0.01 – 0.13 0.034 0.10 0.05 -0.00 – 0.20 0.049
N 0.01 0.04 -0.07 – 0.10 0.732
K 0.07 0.05 -0.03 – 0.17 0.161
CEC -0.09 0.05 -0.20 – 0.01 0.078
Observations 221 221
R2 Nagelkerke 0.056 0.080
Diversity ~ nutrient composition
## Random effect variances not available. Returned R2 does not account for random effects.
[2020]glmmTMB(diversity ~ CEC + K + N, data =dat_t0, ziformula=~1, family = gaussian)
  Best model
Predictors Estimates std. Error CI p
Count Model
(Intercept) 0.71 0.03 0.65 – 0.78 <0.001
CEC -0.07 0.05 -0.17 – 0.02 0.110
K 0.15 0.05 0.05 – 0.24 0.003
N 0.13 0.04 0.05 – 0.21 0.001
(Intercept) 0.50 0.46 – 0.55
Zero-Inflated Model
(Intercept) -21.42 8108.10 -15913.01 – 15870.17 0.998
Observations 223
R2 conditional / R2 marginal NA / 0.015 
## Random effect variances not available. Returned R2 does not account for random effects.
[2021]glmmTMB(diversity ~ CEC + P, data =dat_t1, ziformula=~1, family = gaussian)
  Best model
Predictors Estimates std. Error CI p
Count Model
(Intercept) 1.06 0.04 0.99 – 1.14 <0.001
CEC -0.07 0.04 -0.15 – 0.01 0.068
P 0.11 0.04 0.03 – 0.19 0.006
(Intercept) 0.32 0.27 – 0.38
Zero-Inflated Model
(Intercept) -2.30 0.38 -3.04 – -1.56 <0.001
Observations 86
R2 conditional / R2 marginal NA / 0.013 
## Random effect variances not available. Returned R2 does not account for random effects.
[2022]glmmTMB(diversity ~ N + P + K, data =dat_t2, ziformula=~1, family = gaussian)
  Best model
Predictors Estimates std. Error CI p
Count Model
(Intercept) 0.99 0.05 0.89 – 1.09 <0.001
N 0.14 0.06 0.02 – 0.26 0.026
P -0.14 0.06 -0.25 – -0.04 0.009
K 0.18 0.06 0.06 – 0.30 0.004
(Intercept) 0.45 0.38 – 0.52
Zero-Inflated Model
(Intercept) -22.23 12526.09 -24572.90 – 24528.45 0.999
Observations 84
R2 conditional / R2 marginal NA / 0.028
Composition ~ nutrient composition

2020

## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = ass.rel.t0 ~ N + P + K + C + pH + CEC, data = blocksum[, c(90:95)])
##          Df Variance      F Pr(>F)
## N         1 0.029990 0.9346  0.496
## P         1 0.031152 0.9708  0.454
## K         1 0.029570 0.9215  0.526
## C         1 0.041718 1.3001  0.220
## pH        1 0.031420 0.9792  0.460
## CEC       1 0.037770 1.1771  0.278
## Residual  7 0.224617

## 
## ***VECTORS
## 
##        NMDS1    NMDS2     r2 Pr(>r)
## N   -0.69511 -0.71890 0.1044  0.533
## P    1.00000 -0.00196 0.2460  0.235
## K    0.36640  0.93046 0.1748  0.342
## C   -0.83025  0.55739 0.1548  0.440
## pH   0.88204  0.47117 0.3441  0.103
## CEC  0.11230  0.99367 0.1267  0.522
## Permutation: free
## Number of permutations: 999

2021

## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = ass.rel.t1 ~ N + P + K + C + pH + CEC, data = blocksum[, c(90:95)])
##          Df Variance      F Pr(>F)  
## N         1 0.048054 1.2008  0.268  
## P         1 0.034964 0.8737  0.642  
## K         1 0.027678 0.6917  0.879  
## C         1 0.029522 0.7377  0.813  
## pH        1 0.056124 1.4025  0.096 .
## CEC       1 0.056163 1.4035  0.121  
## Residual  7 0.280120                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## ***VECTORS
## 
##        NMDS1    NMDS2     r2 Pr(>r)  
## N    0.44544  0.89531 0.3457  0.092 .
## P   -0.86862 -0.49548 0.1385  0.406  
## K   -0.75123 -0.66004 0.1268  0.483  
## C    0.84944  0.52768 0.0169  0.922  
## pH  -0.92571  0.37822 0.3009  0.132  
## CEC -0.88035  0.47432 0.0770  0.658  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

2022

## Permutation test for rda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: rda(formula = ass.rel.t2 ~ N + P + K + C + pH + CEC, data = blocksum[, c(90:95)])
##          Df Variance      F Pr(>F)
## N         1 0.056761 1.3466  0.187
## P         1 0.031928 0.7574  0.740
## K         1 0.041588 0.9866  0.470
## C         1 0.032309 0.7665  0.760
## pH        1 0.032999 0.7829  0.724
## CEC       1 0.043354 1.0285  0.459
## Residual  7 0.295060

## 
## ***VECTORS
## 
##        NMDS1    NMDS2     r2 Pr(>r)  
## N   -0.64201 -0.76669 0.0366  0.820  
## P    0.95880  0.28409 0.2378  0.211  
## K    0.61398 -0.78932 0.0906  0.593  
## C   -0.02123 -0.99977 0.3651  0.070 .
## pH   0.94767  0.31924 0.2139  0.278  
## CEC  0.33525 -0.94213 0.1705  0.331  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 999

2020-2022

Interpretation

Q3 Are/how are plant species performances affected by proximity to fallen logs?

Overview of results

2021

There are three response variables: presence (binomial), count (truncated poisson), and per capita biomass (gaussian).

Presence

  • The best fit model for presence is a tie between between physical barrier and the additive model of physical barrier + nutrient island (deltaAIC=0.32). I present the model analysis of the additive model.
  • Both TRCY and TROR are significantly more likely to be present where there is a physical barrier as compared to in the open (TRCY: p=0.03, TROR: p=0.01)
  • Nutrient island does not significantly explain variation in presence/absence in any species. However, there is a trend across all species that they are more likely present in log-legacy plots.

Count

  • The best fit model for count is nutrient island.
  • Nutrient island does not strongly significantly explain variation in presence/absence in any species. However, TRCY is marginally significantly affected by nutrient island, where there are more TRCY in plots where there is a `log legacy’ (there had been a decomposing log in the plot before the experiment began)

Biomass

  • The best fit model for per capita biomass is the nutrient island, then physical barrier island, then the additive model. They are all within 2 AICc points of each other (deltaAIC=1.55 and 1.76, respectively). I present the model analysis of the additive model.
  • Nutrient island significantly explains variation GORO biomass, where there are larger GORO plants in plots where there is a `log legacy’ (p=0.03).
  • Physical barrier does not explain variation in biomass for any species.

2022

There are again three response variables: presence (binomial), count (truncated poisson), and per capita biomass (gaussian).

Presence

  1. The best fit model for presence is a tie between between physical barrier and the additive model of physical barrier + nutrient island (deltaAIC=0.24). I present the model analysis of the additive model.

  2. Both TRCY and TROR are significantly more likely to be present where there is a physical barrier as compared to in the open (TRCY: p=0.007, TROR: p=0.03)

  3. Nutrient island significantly explains the presence/absence of TRCY, where TRCY is more likely to be present in log legacy plots as compared to open legacy plots (p=0.02)

Count

  1. The best fit model for count is a tie between nutrient island and physical barrier. I present the model analysis of the nutrient island model.

  2. Nutrient island does not strongly significantly explain variation in presence/absence in any species.

Biomass

  1. The best fit model for per capita biomass is the nutrient island, then the additive model, then the physical barrier model. They are all within 1 AICc points of each other (deltaAIC=0.31 and 0.46, respectively). I present the model analysis of the additive model.

  2. Nutrient island significantly explains variation in TROR biomass, where log legacy plots have larger TROR as compared to open legacy plots (p=0.04)

  3. Physical barrier significantly explains variation in TROR biomass, where plots in the open have larger TROR than plots near a physical barrier (either log or PVC pipe) (p=0.02)

Statistical Methods

The basic approach is to analyse count and biomass data from 2021 and 2022 sowing experiment. We sowed 15 seeds for each species into 16 plots in each of the seven blocks. There are 6 plot type treatments. The gap and open treatments each have four replicates per block, and the insitu_log, insitu_pvc, open_with_log, and open_with_pvc each have two replicates per block. note: For each treatment of plot type, there is also a dispersal treatment, but I do not analyse that here (yet).
- For count data I use a glmmTMB to run a generalized linear mixed effects model approach to analyse the data. The hurdle model approach I am using is to first code presence/absence as 0 or 1 (1 being any nonzero count value) and run this analysis as a binomial regression. Depending on the model fit and residual dispersion (using DHARMa), I then either run a truncated poisson or truncated negative binomial regression on the count data. Though I did run alternative versions of hurdle models that predict count while accounting for zero inflation, I felt that analyzing the presence/absence, and then analyzing count, may be revealing more of the biology of the system, where presence/absence was distinctly affected by treatments while count was less so.
- For biomass data I use linear mixed effects model to analyse the per capita biomass data. I first do a log(1+n) transformation on the per capita biomass, then analyze. I chose to do per capita biomass because of recalcitrant (!!) residual dispersion in total biomass, even after attemtps at log and square root transformations.
I chose to do a model selection approach, which I’m trying to move away from generally. However, the models I use to analyze the data represent different hypotheses about why plant species might perform differently, and so I chose to do a model comparison and selection approach as a way to not only understand how plant species performance is affected by proximity to fallen logs, but also gain insight as to why this may be the case.

N.B.:Model comparison might not be the approach we want, and we may want to run the analysis differently. I’m open to discussion and change on this point.

After model comparison, I use estimated marginal means (package emmeans) on the best fit model and compute the significance of the difference in estimates.

The models I use test the hypotheses from above: Log decomposition creates islands of fertility directly around the fallen log and fallen logs alter the microclimate directly around them by providing shade.

Analysis

2021
ggarrange(pl1, pl2, pl3, ncol=1, common.legend = T)

2022
ggarrange(pl4, pl5, pl6, ncol=1, common.legend = T)